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1.  Introduction 

To  perform  mesoscale  atmospheric  modeling  on  a  computer,  one  immediately  runs  into  the  problem  of  defining  the  com¬ 
putational  domain.  At  some  point,  there  has  to  be  an  edge  to  the  computational  domain,  but  the  physical  atmosphere  lacks 
any  edges.  How,  then,  can  we  define  a  computational  boundary  where  no  physical  boundary  exists?  The  answer  of  course  is 
to  define  a  non-reflecting  boundary  condition  (NRBC).  How  best  to  define  such  a  boundary  has  been  an  active  area  of  research 
for  approximately  30  years.  Ideally,  an  NRBC  will  be  stable,  accurate,  fast,  and  easy  to  implement;  realistically,  one  must  gen¬ 
erally  choose  two  or  three  of  those  criteria,  at  best. 

There  are  typically  two  approaches  to  NRBC  development.  The  first  is  to  prescribe  the  behavior  at  the  boundaries  in  such  a 
way  as  to  reduce  any  spurious  reflections.  Early  examples  include  the  Sommerfeld-condition-based  work  of  Orlanski  [26] 
and  the  Pade  approximations  of  Engquist  and  Majda  [4,5].  This  approach  was  expanded  by  Higdon  [13-19]  and  subsequently 
automated  by  Givoli,  Neta,  and  van  Joolen  [7-10,29-32].  The  Orlanski  scheme  and  the  Engquist-Majda  scheme  are  less  accu¬ 
rate  than  their  successors;  however,  the  Higdon  scheme  and  its  offshoots  suffer  from  very  high  computational  overhead. 

The  second  approach  is  to  surround  the  domain  with  a  more  dispersive  computational  medium,  so  that  incoming  waves 
enter  the  absorbing  layer  and  diffuse  to  zero  before  their  reflections  re-enter  the  original  domain.  Examples  include  the  per¬ 
fectly  matched  layer  (PML)  developed  by  Berenger  [1  ],  applied  to  the  linearized  shallow  water  equations  by  Navon  et.  al.  [24] 
and  to  the  linearized  Euler  equations  by  Hu  [20-22],  and  the  sponge  layer  used  by  Giraldo  and  Restelli  [6[.  This  approach 
requires  additional  storage  and  computation  time  for  the  expanded  domain,  and  some  reflections  are  still  evident  when 
the  theoretically-exact  absorbing  layer  is  applied  to  a  discrete  computational  domain. 

Here  we  apply  the  Higdon  scheme  to  the  linearized  Euler  equations.  We  take  advantage  of  the  Givoli-Neta-van  Joo¬ 
len  automation  and  make  subsequent  improvements  to  reduce  the  computational  overhead.  This  method  removes 
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approximately  97%  of  the  Sommerfeld-condition’s  reflection  error  with  only  a  modest  increase  to  the  computational 
time. 

The  rest  of  the  paper  is  organized  as  follows:  In  Section  2  we  outline  the  problem  under  consideration,  the  linearized  Euler 
equations  in  2-D  with  no  advection,  solved  in  an  infinite  domain  with  NRBCs  on  all  four  sides.  Section  3  details  the  NRBCs 
and  their  application  to  the  linearized  Euler  equations.  In  Section  4  we  derive  the  Klein-Gordon  equation  from  the  linearized 
Euler  equations  with  no  mean  flow.  We  discuss  the  finite  difference  discretization  for  the  NRBCs  and  the  interior  scheme  in 
Sections  5  and  6,  and  we  provide  a  numerical  example  in  Section  7.  We  then  list  some  areas  for  further  research  (Section  8) 
and  summarize  our  results  (Section  9). 


2.  Problem  statement 


Consider  the  linearized  Euler  equations  in  an  open  domain.  For  simplicity  we  assume  that  the  domain  has  a  flat  bottom 
and  that  there  is  no  advection,  although  this  assumption  may  be  removed  in  future  studies.  A  Cartesian  coordinate  system 
(x,y)  is  introduced,  as  shown  in  Fig.  1. 

The  non-linear  Euler  equations  are 


+  3x(pu)  +  Sy(pv)  =  0, 

dtu  +  udxu  +  vdyii  +  ^dxp=  Jv, 

0tv  +  udxv  +  v6yv  +  yp=  - fu , 

0fP  +  udxp  +  vdyp  +  yp(8xu  +  8yv)  =  0, 


(1) 


where  we  use  the  following  shorthand  for  partial  derivatives 

0  02 
da~8a"  dab~8o8b ’ 

and  t  denotes  the  time,  u(x,y,  t)  and  v(x,y,  t)  the  unknown  velocities  in  the  x  and  y  directions,  p(x,y,  t)  the  density,  p(x,y,  t) 
the  pressure,/  the  constant  Coriolis  acceleration  due  to  the  Earth’s  rotation,  and  y  =  cp/cv  the  constant  ratio  of  specific  heats. 
Linearizing  these  equations  about  mean  zero  velocities,  constant  mean  density  p0  and  constant  mean  pressure  p0  (see,  e.g., 
[20]  or  [23]),  we  get 


0tp  +  po(0xu  +  0yV)  =  0,  (2a) 

8tu  +  ^-8xp=fv,  (2b) 

P  o 

8tv  +  ^-8yP  =  ~fu,  (2c) 

P  o 

8tp  +  yp0(8xu  +  8yv)  =  0.  (2d) 


At  x  — >  oo  the  solution  is  known  to  be  bounded  and  not  to  include  any  incoming  waves.  To  complete  the  statement  of  the 
problem,  initial  values  for  u,  v,  p  and  p  are  given  at  time  t  =  0  in  the  entire  domain. 


Fig.  1.  An  open  domain  Q  truncated  by  artificial  boundaries  rN,  rw,  rs,  and  rE. 
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We  now  truncate  the  infinite  domain  by  introducing  an  artificial  boundary  r,  with  rN  located  at  y  =  yN,  rw  located  at 
x  =  xw,  r s  located  at  y  =  ys,  and  rE  located  at  x  =  xE  (see  dotted  lines  in  Fig.  1 ).  To  obtain  a  well-posed  problem  in  the  finite 
domain  Q  we  need,  instead  of  the  condition  at  infinity,  a  single  boundary  condition  on  each  of  the  artificial  boundaries 
I'n  w.s.e-  This  should  be  a  non-reflecting  boundary  condition  (NRBC).  We  shall  apply  a  high-order  NRBC  for  the  variables, 
as  described  in  the  following  section. 


3.  Higdon-type  NRBCs 

On  the  artificial  boundaries  r  we  use  one  of  the  Higdon  NRBCs  [19].  These  NRBCs  were  presented  and  analyzed  in  a  se¬ 
quence  of  papers  [13,15-18]  for  non-dispersive  acoustic  and  elastic  waves,  and  were  extended  in  [19]  for  dispersive  waves. 
Their  main  advantages  are  as  follows: 


1.  The  Higdon  NRBCs  are  very  general,  namely  they  apply  to  a  variety  of  wave  problems,  in  one,  two,  and  three  dimensions 
and  in  various  configurations. 

2.  They  form  a  sequence  of  NRBCs  of  increasing  order.  This  enables  one,  in  principle  (leaving  implementational  issues  aside 
for  the  moment),  to  obtain  solutions  with  unlimited  accuracy. 

3.  The  Higdon  NRBCs  can  be  used,  without  any  difficulty,  for  dispersive  wave  problems  and  for  problems  in  stratified  media. 
Most  other  available  NRBCs  are  either  designed  for  non-dispersive  media  (as  in  acoustics  and  electromagnetics)  or  are  of 
low  order  (as  in  meteorology  and  oceanography). 

The  scheme  used  here  is  different  than  the  original  Higdon  scheme  [19]  in  the  following  ways: 


1 .  The  discrete  Higdon  conditions  were  developed  in  the  literature  up  to  third  order  only,  because  of  their  algebraic  complex¬ 
ity  which  increases  rapidly  with  the  order.  Givoli  and  Neta  [8]  showed  how  to  easily  implement  these  conditions  to  an  arbi¬ 
trarily  high  order.  The  scheme  is  coded  once  and  for  all  for  any  order;  the  order  of  the  scheme  is  simply  an  input  parameter. 

2.  The  original  Higdon  conditions  were  applied  to  the  Klein-Gordon  linear  wave  equation  and  to  the  elastic  equations.  Here 
we  show  how  to  apply  them  to  the  linearized  Euler  equations  (2b). 

3.  The  Higdon  NRBCs  involve  some  parameters  which  must  be  chosen.  Higdon  [19]  discusses  some  general  guidelines  for 
their  manual  a  priori  choice  by  the  user.  Neta  et.  al.  [25]  showed  how  a  simple  choice  for  these  parameters  can  dramatically 
simplify  the  calculations  and  enable  implementation  of  NRBCs  of  much  higher  order  with  less  computational  overhead. 


The  Higdon  NRBC  of  order  J  is 


Hj  : 


'  1 

+  CjSx) 


j=l 


r\  =  0  on  rE, 


(3) 


where  ;/  represents  any  one  of  the  state  variables  p,  u,  v.  p.  Here,  the  C )  are  parameters  which  have  to  be  chosen  and  which 
signify  phase  speeds  in  the  x-direction.  The  boundary  condition  (3)  is  exact  for  all  waves  that  propagate  with  an  x-direction 
phase  speed  equal  to  any  of  Ci . . .  Cj.  This  is  easy  to  see  from  the  reflection  coefficient  (see  below).  For  the  boundary  rw  we 
replace  0X  by  -Sx.  Likewise,  on  rNS  we  use  ±9y.  Givoli  and  Neta  [8]  and  Dea  et.  al.  [2]  summarize  several  observations  about 
these  NRBCs,  most  of  which  we  omit  here  for  brevity.  One  point  is  worth  repeating:  For  the  Klein-Gordon  equation,  Higdon 
showed  [19]  that  the  reflection  coefficient  for  an  NRBC  of  order  J  can  be  given  by 


i*i=ri 

j=i 


Cj  -  Cx 
Cj  +  cx 


(4) 


where  Cx  is  the  wave  speed  in  the  x  direction.  Note  that  this  reflection  coefficient  is  a  product  of  J  factors,  each  of  which  is  less 
than  one.  Consequently,  increasing  the  order  J  of  the  NRBC  will  result  in  reduced  reflections,  even  if  the  Cj  values  are  sub- 
optimal.  We  will  take  advantage  of  this  fact  when  we  discuss  the  NRBC  discretization  in  Section  5,  where  we  show  that  a 
certain  choice  of  C,  can  reduce  the  computational  complexity  from  0( 3J)  to  0(J2). 


4.  Equivalence  of  linearized  Euler  equations  and  Klein-Gordon  equation 

Higdon  showed  in  [19]  that  this  NRBC  formulation  is  compatible  with  the  Klein-Gordon  (dispersive  wave)  equation 

tfri  -  C20\/2rj  +f2ri  =  0.  (5) 

Hence,  if  we  can  show  that  (2b)  is  equivalent  to  (5),  we  can  claim  that  this  NRBC  formulation  will  be  stable  here. 
Differentiate  (2bd)  with  respect  to  t 


9»P  +  l’Po(exrU  +  6ytV)  =  0. 


(6) 
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Now  differentiate  (2bb)  with  respect  to  x  and  (2bc)  with  respect  to  y  and  add 

dxftf  +  SfyV  +  — —  (pxxP  d-  ^yyP)  =  f  (^x^  ~  dyll). 

P  o 

Now  substitute  (7)  into  (6) 

6  ttP  -  yP[°  (dxxP  +  s  yyP)  +fyp0(dxv  -  dyU)  =  0. 

P  0 

Differentiate  (2bb)  with  respect  to  y  and  (2bc)  with  respect  to  x  and  subtract 

dytu-Qxtv=f(dyv  +  dxu). 


(7) 

(8) 

(9) 


If  we  substitute  this  result  into  the  time  derivative  of  (8),  removing  the  V  x  u  term,  and  apply  (2bd),  we  get  an  equation 
for  8 tp,  which  we  can  integrate  over  time  to  get 

0«p-^V2p+/2(p-po)  =  O,  (10) 

P  o 

which  gives  us  the  Klein-Gordon  equation  for  the  pressure  perturbation  p  -  p0  with  wave  speed  y/yp0/p0. 

5.  Discretization  of  NRBCs 

The  Higdon  condition  Hj  is  a  product  of  J  operators  of  the  form  8f  +  Cfix.  Consider  the  following  finite  difference  approx¬ 
imations  (see  e.g.  [27]): 

i-s - 


0t  — 


i-s ; 

st 


8x 


(11) 


In  (1 1 ),  St  and  5x  are,  respectively,  the  time-step  size  and  grid  spacing  in  the  x  direction,  I  is  the  identity  operator,  and  St 
and  S~  are  backward  shift  operators  defined  by 

«  =  S-Xri;q  =  r,np_hq.  (12) 

Here  and  elsewhere,  tfL  is  the  FD  approximation  of  rj(x,y,  t)  at  grid  point  (xp,y  )  and  at  time  t„.  We  use  (11)  in  (3)  to 


obtain 


n(4f+c- 

u=l  v 


i-s; 

5x 


nh  =  o. 


(13) 


Here,  the  index  E  corresponds  to  a  grid  point  on  the  boundary  rE.  On  the  other  open  boundaries,  the  normal  derivatives 
and  shift  operators  should  be  adjusted  accordingly. 

Givoli  and  Neta  [8]  showed  how  to  implement  the  Higdon  NRBCs  to  any  order  using  a  simple  algorithm.  Their  algorithm 
requires  the  summation  of  0( 3J)  terms.  However,  if  we  make  the  simplification 

Q  =  C0  Vj  €  1 ..  J,  (14) 

then  we  can  simplify  the  summation  to 


J  J-P 


Z=  (al  +  bs;  +  cs;)Jr,«q  =  =  0, 

/M)  y=0 

a  =  1  -  c, 

b  =  - 1, 

st 


where 


(15) 


C  =  -Cn 


a=J-P-y- 

This  summation  consists  of  only  0(J2)  terms,  reducing  the  computational  time  considerably.  As  shown  in  the  discussion  of 
the  reflection  coefficient  in  Section  3,  this  choice  of  C,  will  be  less  accurate  than  an  optimized  selection  of  C/s  based  on  dis¬ 
persive  wave  speeds;  however,  the  results  will  still  improve  as  we  increase  our  order  J.  We  accept  the  slight  loss  of  accuracy 
in  exchange  for  the  significant  increase  in  speed.  Note  also  that  the  J  =  1  case,  with  this  choice  of  Q,  is  the  classic  Sommerfeld 
radiation  condition.  We  choose  this  particular  NRBC  wave  speed  as  a  middle  ground  between  the  phase  speeds  of  the  dis¬ 
persive  waves  and  the  lowered  normal  velocities  associated  with  non-zero  incidence  angles. 


6.  Discretization  in  the  interior 


We  consider  explicit  FD  interior  discretization  schemes  for  the  linearized  Euler  equations  (2b)  to  be  used  in  conjunction 
with  the  Hj  condition.  The  interaction  between  the  Hj  condition  and  the  interior  scheme  is  a  source  of  concern,  since  simple 
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choices  for  an  explicit  interior  scheme  turn  out  to  give  rise  to  instabilities.  The  effort  to  contrive  a  compatible  discretization 
scheme  was  described  in  [2]  for  the  linearized  Euler  equations  without  Coriolis.  There,  we  used  a  one-sided  differencing 
scheme  for  the  interior,  such  that  the  discretized  system  was  equivalent  to  the  standard  second-order  centered-difference 
scheme  for  the  scalar  wave  equation  inp,  which  Higdon  proved  in  [19]  was  compatible  with  the  NRBC  formulation.  However, 
subsequent  work  has  shown  that  adding  the  Coriolis  terms  to  this  scheme  results  in  a  system  which  cannot  be  converted  to 
the  Klein-Gordon  equation.  Hence,  another  approach  is  needed. 

Let  us  reconsider  a  second-order  centered-difference  (leap-frog)  scheme, 


rj(a  +  da)  -  rj(a  -  da) 
2da 


(16) 


where  rj  denotes  any  of  our  four  state  variables,  and  a  denotes  any  of  our  spatial  or  temporal  variables.  Using  the  shift  oper¬ 
ator  notation  from  the  preceding  section,  we  define  our  difference  approximations  as 


ae  {x,y,  t}. 


(17) 


From  this  definition,  we  propose  the  following  discretization  scheme  for  (2b): 


Atp  +  p0(Axu  +  Ayv)  =0, 

(18a) 

AtU  +  -J-Axp  =fv. 

Po 

(18b) 

Atv  +  —  Ayp  =  —fu, 

Po 

(18c) 

Atp  +  yp0(Axu  +  AyV)  =  0. 

(18d) 

Apply  Ax  to  (18bb),  Ay  to  (18bc),  At  to  (18bd),  and  make  the  appropriate  substitution.  This  gives  us 


At A tp  =  ^-{AxAxp  +  AyAyp)  -fyp0(Axv  -  A yu). 
P  0 


(19) 


If/  =  0,  then  this  discretization  is  equivalent  to  a  scalar  wave  discretization.  Hence,  in  the  absence  of  Coriolis  forces,  the  dis¬ 
cretization  scheme  (18b)  is  compatible  with  the  discrete  Higdon  NRBCs  (as  modified  below). 

Continuing  our  derivation,  we  apply  Ay  to  (18bb)  and  Ax  to  (18bc),  then  subtract  and  combine  terms  to  get 


At(Ayu-  Axv)  =f(Ayv  +  Axu). 

We  then  substitute  (18bd)  into  this  result  to  get 

f 

AAAvu  -  Axv)  =  A„. 

YPo 

If  we  apply  At  to  (19)  and  incorporate  (21),  we  get 

YPo, 


At 


At  A  tp  -  (AXA  xp  +  AyAyP)  +f2p 

Po 


=  0. 


(20) 

(21) 

(22) 


Thus,  the  quantity  inside  the  brackets  is  constant  from  one  time  step  to  the  next.  Since  this  equation  applies  to  our  initial 
state,  then  the  quantity  within  the  brackets  must  initially  be  zero  and  thus  remain  zero  always;  hence 


AcAtp  =  ^P-{AxAxp  +  AyAyp)  -f2p. 
Po 


(23) 


If  we  expand  our  A|XJ,  t|  symbols  into  their  corresponding  shift  operators  and  apply  them  to  the  state  variable  p,  we  see 
that  (23)  is  actually 


Pu2  -  2Pij  +  Pu2 


(2dty 


.  YPo  ( Ply  -  2Pij  +  Ply  +  Pij+2  -  2 Pij  +  Pij-2 


Po 


(2dx) 


(2<5y) 


~f  Pi- 


(24) 


This  equation  is  the  standard  second-order  centered-difference  scheme  for  the  Klein-Gordon  equation  on  a  double-sized 
grid.  Hence,  the  appropriate  discretization  for  the  Higdon  scheme  is  not  (11)  but 


—2 


3t  — 


I~Sj 
25 1 


I-S 


-2 


25x 


(25) 


Remark  1.  This  scheme  cannot  resolve  the  shortest  wavelengths  resolvable  by  the  interior  scheme.  Rather,  it  treats  a 
particular  wave  as  an  overlap  of  two  waves,  each  with  twice  the  wavelength  of  the  original  wave,  and  it  resolves  the  waves 
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Fig.  2.  A  computational  wave  (red)  interpreted  as  an  overlap  of  two  longer-wavelength  waves  (green  and  blue).  (For  interpretation  of  the  references  to 
colour  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


on  alternating  time  steps  (see  Fig.  2).  Despite  this  apparent  deficiency,  we  have  observed  no  loss  of  stability  even  in  a 
contrived  example  featuring  such  short  wavelengths. 


7.  Numerical  example 


Let  us  consider  a  simple  numerical  example.  We  look  at  a  square  domain  10  km  on  each  side,  subdividing  it  into  a 
100  x  100  computational  domain  with  the  Higdon-like  NRBCs  on  all  four  sides  (see  Fig.  1).  We  define  the  NRBC  for  all  four 
state  variables,  because  we  use  a  centered-difference  discretization.  A  one-sided  discretization  scheme  would  enable  us  to 
use  a  natural  boundary  condition  on  one  side  or  another,  but  previous  experiments  have  shown  that  such  a  scheme  is 
unstable. 

Using  a  mean  atmoshperic  density  of  1.2  kg/m3  and  pressure  of  1.01xl05N/m2  [12],  a  Coriolis  value  of 
f  =  7.292116  x  10'5  rad/s  [28],  and  zero  advection,  our  initial  condition  is  a  cosine  bubble  in  the  center  of  the  domain 


P°y 

ply 


Po  (  1  + 


:  d  <  r 


Po 


:  otherwise 


( 


Po 


.  d  <  r 
:  otherwise, 


(26) 


where 


d=  \J  (x-xc)2  +  (y-yc)2, 

r  =  1  km, 


and  xc  and  yc  denote  the  center  of  the  domain.  The  initial  condition  for  p  is  chosen  to  maintain  constant  potential  temper¬ 
ature  with  the  pressure  perturbation  [3],  For  comparison,  our  reference  solution  domain  is  30  km  wide  and  30  km  high,  with 
the  domain  of  interest  in  the  center.  We  define  the  normalized  error  norm  for  each  state  variable  rj  as 

£„  =  - - - - - . . . ,  (27) 

^E^EjlV/o(U)2 

where  Nx,  Ny  are  the  number  of  grid  points  in  the  x  and  y  directions,  respectively,  rij  is  a  solution  state  variable  using  the  J- 
order  NRBC,  and  ;/0  is  the  reference  solution.  We  divide  by  the  norm  of  ij0  so  that  the  error  norms  of  each  state  variable  are 
approximately  the  same  order  of  magnitude.  Our  time  step  is  set  to  90%  of  the  maximum  St  allowed  by  the  CFL  limit,  where 
the  CFL  limit  is 


Co 


+  |c«f) 


(28) 


where  C0  =  y/yPo/ p0  is  the  acoustic  wave  speed.  Using  the  discretization  scheme  (18b),  we  run  the  simulation  up  to  t  =  24, 
long  enough  for  the  primary  wave  to  exit  the  computational  domain  with  the  wave  trough  just  passing  through  the  corners. 
Figs.  3-6  show  the  four  state  variables  at  the  end  of  the  run  for  J  =  10.  Fig.  7  compares  the  state  variable  u  at  the  end  of  the 
J  =  1  and 7=10  cases;  note  the  dramatic  reduction  in  spurious  reflection  in  the 7=10  case.  (We  have  observed  here  that  if 
the  order  is  increased  beyond  7  =  10,  numerical  instabilities  appear  which  destroy  the  solution.  We  observed  the  same  insta¬ 
bility  for 7  >  10  in  developing  the  numerical  examples  of  [25]  for  the  Klein-Gordon  equation.  The  cause  is  unknown  but  is 
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Density  perturbation  x  ^g4 
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Density  perturbation,  reference  solution 


50  100 

Error  norm-0.029239 
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Fig.  3.  The  solution  for  the  density  p  using  ]  =  10.  (TL)  Computed  solution  using  NRBCs.  (Right)  Reference  solution.  (CL)  Reference  solution  domain 
corresponding  to  NRBC  solution  domain.  (BL)  Delta  between  computed  and  reference  solutions,  with  error  norm  computed  by  (27). 


X-direction  velocity  perturbation 
100 


X-direction  velocity  perturbation,  reference  solution 


Truncated  reference  solution 

100  F 


50  100 

Error  norm-0.062476 


200 


Fig.  4.  The  solution  for  the  x-direction  velocity  u  using  J  =10.  (TL)  Computed  solution  using  NRBCs.  (Right)  Reference  solution.  (CL)  Reference  solution 
domain  corresponding  to  NRBC  solution  domain.  (BL)  Delta  between  computed  and  reference  solutions,  with  error  norm  computed  by  (27). 


suspected  to  be  due  to  round-off  errors  associated  with  small  individual  values  in  the  summation  of  the  high-order  deriv¬ 
atives.)  Table  1  and  Fig.  8  shows  the  error  norms  (27)  for  each  state  variable  as  J  goes  from  1  to  10. 
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Y-direction  velocity  perturbation 
100 


200 
0.05 


0 

150 

-0.05 


50  100 


Truncated  reference  solution 
100 


100 
0.05 

0  50 

-0.05 


50  100 

Solution  delta 


x  10 
4 


-50 


50  100 


Y-direction  velocity  perturbation,  reference  solution 


-50  0  50  100  150  200 


Error  norm~0. 062477 

Fig.  5.  The  solution  for  the  y-direction  velocity  v  using  J  =  10.  (TL)  Computed  solution  using  NRBCs.  (Right)  Reference  solution.  (CL)  Reference  solution 
domain  corresponding  to  NRBC  solution  domain.  (BL)  Delta  between  computed  and  reference  solutions,  with  error  norm  computed  by  (27). 


One  note  about  the  time  step:  we  have  found  experimentally  that  if  St  is  set  to  exactly  the  CFL  limit,  then  the  error 
norm  for  J  =  10  is  only  approximately  55%  the  J  =  1  error  norm.  On  the  other  hand,  if  we  use  exactly  half  the  CFL  limit 
for  our  St,  then  the  error  norm  reduction  is  on  the  order  of  99.5%.  We  chose  90%  of  the  maximum  St  as  a  compromise 
between  increased  accuracy  and  increased  time  step  size.  The  reason  for  the  improved  NRBC  performance  with  smaller 
St  is  unknown,  but  it  may  be  similar  to  the  “safety  factor”  recommended  by  Tannehill  et  al.  for  the  non-linear  Navier- 
Stokes  equations  (see  p.  627  of  [27]).  Analysis  of  this  particular  aspect  of  the  NRBC  performance  will  be  addressed  in 
future  research. 


8.  Areas  for  further  research 

The  preceding  example  demonstrates,  in  a  limited  setting,  that  high-order  Higdon  NRBCs  are  compatible  with  the  line¬ 
arized  Euler  equations.  However,  there  are  far  more  areas  to  explore  in  this  implementation.  The  following  list  shows  some 
of  the  areas  available  for  future  research,  some  of  which  are  currently  under  investigation  by  the  authors: 

1 .  Thorough  investigation  of  the  long-time  stability  for  large  J. 

2.  Extending  the  scheme  to  the  case  of  the  linearized  Euler  equations  with  Coriolis  and  nonzero  mean  flow  (advection). 

3.  Thorough  investigation  of  the  relationship  between  the  time  step  size  (the  Courant  number)  and  the  NRBC  performance 
for  large  J. 

4.  Extending  the  scheme  to  the  full  3-D  system,  including  the  effects  of  gravity. 

5.  Implementing  the  scheme  with  auxiliary  variables,  using  high-order  finite  differences  and  finite  elements,  using  both  the 
Givoli-Neta  AV  formulation  [10]  and  the  Hagstrom-Warburton  variation  [11]. 

6.  Extending  the  scheme  to  permit  incoming  waves,  for  example,  in  a  nested  mesoscale  model. 

7.  Experimenting  with  the  use  of  the  NRBC  with  the  non-linear  Euler  equations  (1)  in  the  computational  domain.  (Need  to 
find  a  stable  interior  scheme-NRBC  combination.) 

9.  Conclusion 

In  this  paper,  we  have  shown  that  Higdon-type  NRBCs  are  compatible  with  the  linearized  Euler  equations  with  Coriolis 
and  zero  mean  flow.  These  NRBCs  provide  greater  accuracy  (reduced  spurious  reflection)  than  the  basic  Sommerfeld  bound¬ 
ary  condition.  A  prototypical  implementation  was  developed,  and  a  numerical  example  demonstrating  the  capabilities  of  the 
scheme  was  provided. 
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Fig.  6.  The  solution  for  the  pressure  p  using  7=10.  (TL)  Computed  solution  using  NRBCs.  (Right)  Reference  solution.  (CL)  Reference  solution  domain 
corresponding  to  NRBC  solution  domain.  (BL)  Delta  between  computed  and  reference  solutions,  with  error  norm  computed  by  (27). 
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Fig.  7.  A  comparison  of  solutions  for  u  using7  =  1  and  7=10.  (TL)  Computed  solution  using7  =  1.  (TR)  Computed  solution  using7  =10.  (BL)  Delta  between 
reference  solution  and  7=1  solution,  with  error  norm  computed  by  (27).  (BR)  Delta  between  reference  solution  and  7  =  10  solution,  with  error  norm 
computed  by  (27). 
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Table  1 

Error  norms  (27)  for  J  e  1 ...  10  with  discretization  scheme  (18b). 


] 

EP 

Eu 

£v 

Ep 

1 

1.5191 

2.0917 

2.0917 

1.5205 

2 

0.42052 

0.61777 

0.61777 

0.42092 

3 

0.18953 

0.30055 

0.30054 

0.18971 

4 

0.11677 

0.19766 

0.19766 

0.11689 

5 

0.081815 

0.14588 

0.14588 

0.081893 

6 

0.061569 

0.11564 

0.11564 

0.061628 

7 

0.048183 

0.095798 

0.095797 

0.04823 

8 

0.03908 

0.082285 

0.082284 

0.039118 

9 

0.033036 

0.071617 

0.071617 

0.033067 

10 

0.029239 

0.062476 

0.062477 

0.029267 

Density  perturbation  error  norm  X-velocity  perturbation  error  norm 


Y-velocity  perturbation  error  norm  Pressure  perturbation  error  norm 


Fig.  8.  Logarithmic  plot  of  state  variable  error  norms  (27)  for  J  e  1 ...  10  with  discretization  scheme  (18b).  (TL)  Error  norms  for  p.  (TR)  Error  norms  for  u. 
(BL)  Error  norms  for  v.  (BR)  Error  norms  forp. 
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